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Abstract. In numerical computations of astrophysical 
flows one frequently has to deal with situations which re¬ 
quire the use of a corotating coordinate system. Then, in 
addition to the centrifugal force, the Coriolis terms have 
to be included in the Navier-Stokes equations. In usual 
numerical treatments of solving the equations of motion 
these non-inertial forces are included explicitly as source 
terms. In this letter we show that this procedure, if ap¬ 
plied to the angular momentum equation, may lead to 
erroneous results. Using the specific sample calculation of 
a protoplanet embedded in a protostellar disk, we then 
demonstrate that only a conservative treatment of the 
Coriolis term in the angular momentum equation yields 
results which are identical to those obtained by computa¬ 
tions performed in the inertial frame. 
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1. Introduction 

In numerical computations of the hydrodynamic evolution 
of astrophysical flows the usage of a rotating coordinate 
system is often desirable. For example, in calculating the 
evolution of an accretion disk in a cataclysmic variable 
system usually a coordinate system which rotates with 
the orbital speed of the secondary is used (eg. Heemskeerk, 
1994). Then the secondary star remains at a fixed location 
in the grid, simplifying the numerical procedure. Also, in 
computations of a disk around a binary star (Artymowicz 
& Lubow, 1994) or with an embedded planet in the disk 
(Kley, 1998, Bryden et al. 1998) corotating coordinates are 
very useful. In particular, this simplifies grid refinement in 
the vicinity of an embedded planet to resolve the detailed 
flow within the planet’s Roche lobe. 

In standard, recent approaches to model the tidal effect 
on an accretion disk in a binary system the Coriolis terms 
are treated as explicit source terms in the equations of 
motion (eg. Heemskerk 1994, Godon 1996, Kunze et al. 
1997). Using this procedure, the Coriolis force does not 
conserve total angular momentum. We are not aware that 


in any of these computations comparison runs between the 
inertial and corotating frame have been performed. As will 
be shown in this letter, using the sample calculation on 
an embedded planet within an accretion disk, the explicit 
approach yields incorrect results concerning the overall 
evolution of the system. We then proceed to demonstrate 
that, on the contrary, the conservative formulation of the 
Coriolis term in the angular momentum equation leads 
to results which are identical to those performed in the 
inertial frame. 

In the following section 2 the physical setup of the 
calculation is described. In section 3, the dependence of 
the obtained results on the treatment of the Coriolis force 
are displayed, and in section 4 we summarize. 

2. Physical Model 

In an accretion disk where the vertical thickness H is usu¬ 
ally much smaller than the distance r from the center 
(H/r « 1) it is customary to work only with vertically 
averaged quantities. Hence, the system under considera¬ 
tion consists of an infinitesimally thin disk with an embed¬ 
ded protoplanet. We work in a rotating coordinate system 
where the planet is at a fixed location. We shall work in 
cylindrical coordinates (r, <p, z), where r is the radial coor¬ 
dinate, ip is the azimuthal angle, and 2 is the vertical axis, 
and the rotation is around the 2 -axis, i.e ft = (0,0, U). 
The origin of the coordinate system is located at the cen¬ 
ter of mass of the system. 

In a flat disc located in the 2 = 0 plane, the velocity 
components are u = (u r ,u v ,0). In the following we will 
use the symbol v = u r for the radial velocity and uj = u^/r 
for the angular velocity of the flow, which are measured 
in the corotating frame. Then the vertically integrated 
equations of motion are 
r)T 

—+ V(£u) = 0, (1) 

^> + V(E™) = Er( u + Si) a -g-E^ + / r (2) 

2^1 + V(EA,u) = -2S„ O r - | - e|£ + /, (3) 
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Fig. 1. 3D-mountain plot of the surface density using a 
rotating coordinate system and the conservative treatment 
(Eq. [|) of the angular momentum equation 


Fig. 2. 3D-mountain plot of the surface density using a 
rotating coordinate system and the explicit treatment (Eq. 
P) of the angular momentum equation 


Here E denotes the surface density and p the vertically 
integrated (two-dimensional) pressure. The gravitational 
potential created by the protostar with mass M s and 
the planet having mass m p is given by 

GM S Gm p 
l r r s | |r — r p | ’ 

where G is the gravitational constant and r s , r p are the ra¬ 
dius vectors to the star and planet, respectively. To model 
the potential of the planet within the Roche lobe we in¬ 
troduce a smoothing length of 1/5 of the Roche radius. 
The effects of viscosity are contained in the terms f r , and 
which give the viscous force per unit area acting in the 
radial and azimuthal direction. The explicit form of these 
terms is given in Kley (1998). 

Note, that in the angular momentum equation (§ the 
Coriolis forces are included explicitly as source terms on 
the right hand side. If treated numerically, they do not 
conserve total angular momentum, and to demonstrate 
the effect of this formulation on the physical results the 
following, additional conservative formulation is used in 
the computations 


<9[E r 2 (uj + fi)] 
dt 


dp 


.< 9 $ 


+ V[Er V + fi)u] = --^-E— + ftp.(4) 


As will be shown, only this latter form leads to correct 
results in the corotating frame. An energy equation is 
not required since an isothermal equation of state is used 
where the pressure is given by p — E . The sound speed 
is set to a constant fraction of the Keplerian velocity 
c s = 0.05uKep, which is equivalent to a fixed relative thick¬ 
ness of the disk, H/r = 0.05. The kinematic viscosity co¬ 
efficient in these test calculations is set to a constant value 



Fig. 3. 3D-mountain plot of the surface density using a 
non-rotating coordinate system (f l = 0) 

v = 10 -5 in dimensionless units (see below), and the bulk 
viscosity is set to zero. For the mass ratio of the planet to 
the star we assume q = m p /M s = 10 -3 . 

2.1. Numerical issues, initial and boundary conditions 

To solve the equations of motion dimensionless units are 
used where the unit of length is given by the distance a 
of the planet to the star. The unit of time is given by 
to = Qp 1 where S7 p is the orbital speed of the planet; 
identical to the speed of the rotating coordinate system. 
The other units follow from these. The disk is non self- 
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Table 1. Numerical parameter of the models 

Model Code Description 

1 Nirvana Rotating frame (conservative, Eq. [Ji 

2 Nirvana Rotating frame (explicit, Eq. |j) 

3 Nirvana Inertial frame (f2 = 0) 

4 WK Inertial frame (12 = 0) 

as 

E 

gravitating, and the density scales out of the problem. ~ 

The time in the result section, and the plots is given units "g 

of the period P p = 27rf2” 1 of the planet. a 

The detailed description of the physical model used > 
for the present analysis is outlined in Kley (1998) and 
is not repeated here. The main results of this letter are 
obtained with a code called Nirvana (Ziegler & Yorke, 

1997), which is a 3D-multi-grid MHD code. It was adapted 
for this 2D, single grid, purely hydrodynamic application 
and extended by including the viscous terms explicitly. 

The code uses a spatially second order accurate, explicit 
method which handles advection by a monotonic trans¬ 
port algorithm. Hence, the advective parts of the code 
conserve globally mass and angular momentum. An ad¬ 
ditional run, using identical physical parameter, has been 
performed with a different, independent code WK devel¬ 
oped by Kley (1989), which uses an implicit treatment of 
viscosity. 

The initial setup consists of a Keplerian disk, with a 
density profile E oc r -1 / 2 which follows in case of a con¬ 
stant kinematic viscosity. Additionally, an initial gap (low¬ 
ering of the density) is assumed around the planet whose 
size can be estimated by equating viscous and gravita¬ 
tional torques (Lin & Papaloizou, 1993). The initial E(r) 
profile is given for example in Fig. 4. The radial velocity 
is set to zero initially. The inner and outer boundary are "g 
closed throughout the evolution. The computational do- « 
main has the extent r m ; n = 0.4, r max = 2.5, and covers > 
the whole azimuthal range in = 0, </J max = 2n. The grid 
resolution consists of 128 x 128 equidistantly spaced grid 
points. The physical parameter for all 4 models presented 
are identical, and we state in Table 1 the main, differing 
computational features of the models. 

3. Numerical results 

This initially axially symmetric disk setup is evolved un¬ 
der the influence of the perturbing potential of the planet, 
which very soon creates non-axisymmetric features in the 
density distribution. Trailing spiral arms are developing 
inside and outside of the planet’s radial location. For de¬ 
tails and the influence of the variation of the physical pa¬ 
rameter see Kley (1998), and Bryden et al. (1998). Here 
we are only interested in numerical aspects of the com¬ 
putations. The total evolution time of the system covers 
about 200 orbital periods of the planet. 

In Figs. 1 to 3 mountain plots of the surface density 
distribution at t=200 are given for the different calcula- 



Fig. 4. The azimuthally averaged surface density for 
model 1. 



Radius 

Fig. 5. The azimuthally averaged surface density for 
model 2. 


tions, where for clearity only half the grid points are dis¬ 
played. While Figs. 1 and 2 refer to the computation using 
a corotating coordinate system Fig. 3 refers to a compu¬ 
tation in the inertial frame where Q has been set to zero. 
In obtaining Fig. 1 the conservative treatment ([|) of the 
angular momentum equation has been used, while Fig. 2 
was computed using the explicit version (i3jj). 
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Clearly, there are distinct differences in the density dis¬ 
tribution in the two formulations of the problem. When 
the explicit form (J3j) of the Coriolis terms is used the ob¬ 
tained density profile is much shallower radially around 
the planet (located at r = 1 and p = n). More matter 
seems to be pushed away from the planet to larger radii. 
This indicates an additional unphysical transport of angu¬ 
lar momentum from inside out in the vicinity of the planet. 
In contrast, the inertial frame calculation (Fig. 3) and the 
first, conservative calculation (Fig. 1) agree very well with 
each other. As the physical results should not depend on 
any assumed rotation of the coordinates we conclude that 
only in the case where the Coriolis terms are treated in a 
conservative fashion the correct results are obtained. 

This important finding is corroborated by looking at 
the time evolution of the radial distribution of the surface 
density averaged over the azimuthal angle <p as given in 
Fig 4 and 5. In the conservative treatment (Fig. 4) the 
averaged density distribution reaches a quasi-stationary 
state which does not evolve any further. This is to be ex¬ 
pected once gravitational and viscous torques are in equi¬ 
librium. In contrast, the explicit formulation shows a con¬ 
tinuing evolution of the density profile in the outer region 
(Fig. 5), where mass appears to be accumulating near the 
outer boundary. 

A physically interesting question is related to possi¬ 
ble mass accretion through the gap to form more massive 
planets as observed in extrasolar systems (Boss, 1996). 
To model such a situation, mass has been removed con¬ 
tinually within the Roche lobe of the planet with a half 
emptying time of about 1/4 orbital periods of the planet 
(Kley 1998). In figure 6 the obtained mass accretion rate 
for the planets (in dimensionless units) is plotted versus 
time for the four models where the model named WK 
has been obtained in the inertial frame with a different 
code (Kley 1989). Obviously the codes Nirvana and WK 
agree very well with each other, as do the conservative 
and inertial frame calculation. Only the explicit form ([j|) 
deviates again strongly and shows additionally a decline 
with time. The reduced accretion onto the planet appears 
to be caused by artificial angular momentum increase of 
the explicit model. 

4. Local analysis 

To verify the conjecture above, we perform a local anal¬ 
ysis of the finite difference equations. Specifically, let us 
consider the change of angular momentum A J within one 
gridcell during one timestep At. For simplicity, axial sym¬ 
metry is assumed, and pressure and viscous effects are 
not considered. The change in the i th radial gridcell with 
volume Vi using the conservative treatment is then given 
purely by the advection term 

Ajcons = _ [ F . +1 ( w + n)a l+ 1 - Fi(u + f?)oj] At, (5) 



Fig. 6. The accretion rate onto the planet versus time for 
the different models. 

with the angular momentum flux (flow per unit length 
and time) Fi{u> + ft) = rf Vi {pJi-i + II) S*_i, through the 
curcumference at = 2irri. To demonstrate the effect, we 
consider a simple upwind scheme on a staggered, equidis¬ 
tant grid, where and u>i are located at Vi + i/ 2 in the 
middle of the cell interfaces r,; and ?y + 1 . The radial veloc¬ 
ities Vi are located at r i: and in (^) positive velocities are 
assumed. For the explcicit treatment (Eq. ^]) we obtain 

AJ? xpl _ _ [F z+1 (u)a i+ i - Fi{u)a,i ] At, 

+ AJ^ + r^QAMi ( 6 ) 

where 

AJf or = -Ej {Vi + u i+1 )Hr l+1/2 Af AK, (7) 

is the change in angular momentum caused by the explicit 
Coriolis term, and AM; is the change of mass in the grid¬ 
cell. Taking now the difference AJ = AJ^ xpl — AJf ons , 
and assuming 

£ = E 0 r s , v = v 0 r q , 

and a quasi-stationary density distribution A Mi = 0, then 
we find 

AJ = 27rf2Af£ 0 no r q+s+3 2x{q + s + 1), (8) 

where x = Ar/?' i+1 / 2 . Thus, for positive velocities which 
inevitably occur at the outer edge of the gap and a posi¬ 
tive gradient of the density (see Fig. 4 and 5), we find that 
AJ > 0. Hence, the explicit treatment shows an artificial 
enhancement of angular momentum, which tends to widen 
the gap created by the planet even further (see Fig.5). We 
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notice, that the identical relation (||) is obtained for neg¬ 
ative velocity (vq < 0) which means that for a standard 
accretion disk with q = — 1 and s negative, the explcicit 
treatment yields to an artificial increase of angular mo¬ 
mentum as well. Both effects play a role in the presented 
model calculations of a planet embedded in an accretion 
disk. 

5. Conclusions 

We have performed numerical computations of a planet 
embedded in a disk using different treatments of the Cori¬ 
olis force in the angular momentum equation. 

Three different models, one with an explicit Coriolis 
force, the second using a conservative formulation for the 
Coriolis term, and the third in the inertial frame (Q = 0) 
have been calculated. It was seen that only the second and 
the third model yield identical results, which were also 
confirmed by an additional calculation using a completely 
independent numerical code. 

We conclude that in performing numerical models 
which use a rotating reference frame, care must be taken 
to incorporate the Coriolis terms correctly in the angular 
momentum equation. As often in numerical computations 
a conservative formulation of the Coriolis force turned out 
to be superior over an explicit, non-conservative treat¬ 
ment. In the future numerical modelers should be aware of 
this potential problem when calculating applications such 
as disks in binary systems, circumbinary disks, rotating 
stars, or similar problems. 
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